%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Titulo: Seguimiento Lagrangiano de Particulas
%
% Autor: Luis Miguel de la Cruz Salas
%
%------------------------ COMPILACION -----------------------------
%latex flowviz.tex
%latex flowviz.tex
% dvips -Ppdf -G0 flowviz.dvi -z -o flowviz.ps
%ps2pdf13 -sPAPERSIZE=letter -dMaxSubsetPct=100 -dCompatibilityLevel=1.3 -dSubsetFonts=true -dEmbedAllFonts=true flowviz.ps
%-----------------------------------------------------------------
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%-----------------------------------------------------------------------
% Letra de 12pt, impresion por los dos lados, tama\~no carta, tipo libro
%-----------------------------------------------------------------------
\documentclass[12pt,twoside,letterpaper,colordvi]{article}

%-----------------------------------------------------------------------
% Packages 
%-----------------------------------------------------------------------
\usepackage[dvips]{graphicx}
\usepackage{epsfig}
\usepackage{spanish}
\usepackage{fancyheadings}
\usepackage{hhline}
\usepackage{moreverb}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{color}
\usepackage[dvips]{hyperref}
%\usepackage{multicol}
%\usepackage{supertab}
\hypersetup{
	pdfauthor={Luis M. de la Cruz Salas},
	pdftitle={Partrack, DCI-DGSCA-UNAM},
	backref,
	breaklinks,
	colorlinks=true
}
%-----------------------------------------------------------------------
% Color Definition
%-----------------------------------------------------------------------
\definecolor{rosa}{rgb}{1,0.5,0.5} 
\definecolor{brown}{rgb}{0.64,0.16,0.16}
\definecolor{ForestGreen}{rgb}{0.13,0.54,0.13}
\definecolor{purple}{rgb}{0.62,0.12,0.94}
\definecolor{DodgerBlue}{rgb}{0.11,0.56,0.98}
\definecolor{RoyalBlue}{rgb}{0.25,0.41,0.88}
\definecolor{gold}{rgb}{0.85,0.66,0}
\definecolor{red}{rgb}{1, 0, 0}
\definecolor{green}{rgb}{0, 1, 0}
\definecolor{capitulo}{rgb}{0, 0, 0.50} % DarkBlue
\definecolor{seccion}{rgb}{0.62,0.12,0.94} % purple
\definecolor{subseccion}{rgb}{0.11,0.56,0.98} %DodgerBlue
%{\textcolor{sectionColor}{ Texto en color } } 

%-----------------------------------------------------------------------
% Text dimensions 
%-----------------------------------------------------------------------
\setlength{\textheight}{22cm}
\addtolength{\textwidth}{56pt} % Aproxx. ~16cm
%
% Becuase I don't need marginal notes, then: 
%
\setlength{\marginparsep}{0pt}
\setlength{\marginparwidth}{0pt}
%
% Header is equal to \textwdith
% This is used when \marginparsep and \marginparwith are not zero
%
%\addtolength{\headwidth}{\marginparsep}
%\addtolength{\headwidth}{\marginparwidth}
%
% Header is equal to \textwdith
%
\setlength{\headwidth}{\textwidth}
%
% We have to substrac some points on even pages,
% because we add some points to \textwidth
%
\addtolength{\evensidemargin}{-56pt}
%\addtolength{\oddsidemargin}{-10pt}

\setlength{\voffset}{-1cm}

%-----------------------------------------------------------------------
% Style of page and headings
%-----------------------------------------------------------------------
\pagestyle{fancyplain}
\renewcommand{\sectionmark}[1]
	{\markright{\thesection\ #1}}
\lhead[\fancyplain{}{\bf\thepage}]
	{\fancyplain{}{\bfseries\rightmark}}
\rhead[\fancyplain{}{\bfseries\leftmark}]
	{\fancyplain{}{\bf\thepage}}
\rfoot{\small\textsf{\copyright 2004-2005 Depvis,DCI-DGSCA}}
\lfoot{\small\textsf{luiggi}}
\cfoot{}

\newcommand{\bmp}[2][t]{\begin{minipage}[#1]{#2}}
\newcommand{\emp}{\end{minipage}}
 

%-----------------------------------------------------------------------
% Title and author
%-----------------------------------------------------------------------
\title{Seguimiento Lagrangiano de Part\'{\i}culas}
\author{Luis M. de la Cruz \\ luiggi@ixtli.unam.mx \\
www.labvis.unam.mx/luiggi}

\begin{document}

\maketitle

\tableofcontents

\section{Descripci\'on y teor\'{\i}a}

Una trayectoria est\'a compuesta por las posiciones subsecuentes que ocupa 
una part\'{\i}cula cuyo movimiento es promovido por un campo vectorial, 
v\'ease figura \ref{figtracer}.

\begin{figure}[h]
\begin{center}
\epsfig{file=figtracer.eps, width=14cm}
\end{center}
\caption{Trayectoria de una part\'{\i}cula inmersa dentro de un campo
vectorial.}
\label{figtracer}
\end{figure}

El seguimiento de part\'{\i}culas se basa en la soluci\'on a la siguiente 
ecuaci\'on diferencial ordinaria:

\begin{equation}\label{ode}
\frac{d \vec{x}(t)}{d t} = \vec{U}(\vec{x}(t), t)
\end{equation}

\noindent
donde $\vec{x}$ representa la posici\'on de una part\'{\i}cula en el instante
de tiempo $t$, y $\vec{U}(\vec{x}, t)$ representa a un campo vectorial,
(velocidad del flujo, vorticidad, fuerza, etc.). Una posici\'on inicial 
$\vec{x}(t_0) = \vec{x}_0$ provee la condici\'on inicial con la que es
posible resolver la ecuaci\'on (\ref{ode}). La soluci\'on de (\ref{ode})
ser\'a una secuencia de posiciones ($\vec{x}_0, \vec{x}_1, \vec{x}_2,
\dots$) que definen una trayectoria. 

El campo vectorial puede provenir de una funci\'on, una simulaci\'on 
num\'erica o una medici\'on experimental.

El dominio espacial se representa mediante una malla que puede ser una 
combinaci\'on de regular o irregular, estructurada o no estructurada, 
uniforme o no uniforme. Ejemplos de diferentes tipos de mallas se muestran
en la figura \ref{figmallas} (v\'ease tambi\'en \cite{}).

\begin{figure}[h]
\begin{center}
\epsfig{file=figmallas.eps, width=15cm}
\end{center}
\caption{Diferentes tipos de mallas: (a) regular, estructurada, uniforme,
(b) regular, estructurada, no uniforme, (c) regular, no estructurada,
(d) irregular, estructurada y (e) irregular, no estructurada.}
\label{figmallas}
\end{figure}

Una vez que se tiene definida la malla y el campo vectorial sobre dicha 
malla, entonces se inicia el c\'alculo de la trayectoria de part\'{\i}culas.
Existen diferentes tipos de trayectorias, las cuales se describen en las 
siguientes secciones.

\subsection{Streamlines}
Las streamlines son curvas que en cada uno de sus puntos
son tangentes al campo vectorial. Esta es una propiedad que proviene de 
la definici\'on de la ecuaci\'on (\ref{ode}).
Las streamlines son \'utiles para analizar campos vectoriales
que no cambian con el tiempo. La figura \ref{stream} muestra ejemplos de
streamlines.

\begin{figure}[h!]
$$\includegraphics[height=8cm]{stream.eps}$$
\caption{Ejemplo de algunas streamlines.}
\label{stream}
\end{figure}

Debido a que en este caso el campo vectorial no var\'{\i}a con el tiempo, la 
ecuaci\'on (\ref{ode}) se puede escribir como

\begin{equation}\label{streamode}
\frac{d \vec{x}}{d s} = \vec{U}(\vec{x}(s)) 
\end{equation}

\noindent
donde $s$ es un par\'ametro que define el inicio y final de la trayectoria.
Si se requiere que la streamline inicie en el punto $\vec{x}_0$, entonces
la condici\'on inicial debe ser $\vec{x} = \vec{x}_0$ en $s = 0$,
en donde, el par\'ametro $s \in R$ y controla la extensi\'on de la
trayectoria. Para resolver la ecuaci\'on (\ref{streamode}), se debe realizar 
una integraci\'on sobre $s$:

\begin{equation}\label{streamintegra}
\int_s^{s+\Delta s} \frac{d \vec{x}}{d s} ds = 
\vec{x}_{s+\Delta s} - \vec{x}_s =
\int_s^{s+\Delta s} \vec{U}(\vec{x}(s)) ds
\end{equation}


Ahora, si suponemos que $\vec{U}$ es constante en el intervalo $[s,
s+\Delta s]$ entonces la integral de la derecha de la ecuaci\'on
(\ref{streamintegra}) se reduce a

\begin{equation}\label{streamintegrada}
\vec{x}_{k+1} = \vec{x}_k + \vec{U}(\vec{x}_k) \Delta s
\end{equation}

\noindent
donde hemos utilizado la notaci\'on $\vec{x}_{s+\Delta s} = \vec{x}_{k+1}$ y
$\vec{x}_{s} = \vec{x}_{k}$. La f\'ormula (\ref{streamintegrada}) se conoce
como la f\'ormula de Euler y es el m\'etodo m\'as simple para integrar la
ecuaci\'on (\ref{streamode}).

Se observa que para construir una streamline, primero se da una posici\'on
inicial en un cierto instante: $\vec{x}(s_0) = \vec{x}_0$. Luego, de
acuerdo con (\ref{streamintegrada}), la posici\'on siguiente en la
trayectoria ser\'a $\vec{x}_1 = \vec{x}_0 + \vec{U}(\vec{x}_0) \Delta s$.
Este proceso se continua tanto como se desee. En general, la velocidad
$\vec{U}(\vec{x}_k)$ no cae en lugares de la malla donde la velocidad es
conocida, por lo tanto, se debe usar una f\'ormula de interpolaci\'on para
calcular la velocidad en cada posici\'on, v\'ease la figura \ref{streams}. 
Cuando el campo vectorial est\'a dado por una funci\'on, no es necesario
realizar ning\'un tipo de interpolaci\'on.

\begin{figure}[h!]
$$\includegraphics[height=8cm]{streamcalc.eps}$$
\caption{C\'alculo de las posiciones de la streamline.}
\label{streams}
\end{figure}

En principio $\Delta s$ debe ser infinitesimal. Sin embargo, en nuestra
simulaci\'on num\'erica estamos limitados por la precisi\'on de la
m\'aquina. Por lo tanto se debe elegir una $\Delta s$ que permita obtener
streamlines que siguan las leyes f\'{\i}sicas y, que por otro lado, sea
posible calcularlas en tiempos relativamente cortos. Una manera de obtener
una $\Delta s$ adecuada es como sigue:

	\begin{enumerate}
	\item Se elige una $\Delta s_1$ inicial.
	\item Se calcula la trayectoria con $\Delta s_1$.
	\item Se elige otra $\Delta s_2$, con la condici\'on 
	$\Delta s_2 < \Delta s_1$.
	\item Se calcula una nueva trayectoria usando $\Delta s_2$.
	\item Se mide la diferencia entre las dos trayectorias. Si la
	diferencia es m\'{\i}nima, dentro del rango de error permitido,
	entonces usamos $\Delta s_1$. En otro caso, renombramos a 
        $\Delta s_2$ como $\Delta s_1$ y repetimos el proceso
	hasta obtener una $\Delta s$ adecuada. 
	\end{enumerate}


El siguiente es un algoritmo general para calcular streamlines:

\strut

%\framebox{ \bmp{5.5in}

\begin{center}\fbox{\bf{Algoritmo 1: Streamlines}}\end{center}

{\small{
\begin{enumerate}
\item Determinar lo siguiente:
	\begin{enumerate}
	\item Tipo de malla sobre la que est\'a definido el campo
	vectorial $\rightarrow$ \textsf{Mesh}
	\item N\'umero de pasos para calcular la trayectoria 
	$\rightarrow$ \textsf{Nsteps}
	\item Tama\~no de la $\Delta s$ para realizar los c\'alculos de la
	trayectoria $\rightarrow$ \textsf{Dstep}
	\item Posici\'on inicial de la part\'{\i}cula 
	$\rightarrow$ \textsf{Position}
	\end{enumerate}
\item Obtener el campo vectorial estacionario $\rightarrow \vec{U}$
\item {\bf\textsf{WHILE}} ( \textsf{Position} $\in$ \textsf{Mesh}
			      \quad {\bf\textsf{OR}} \quad 
			    \textsf{s} $\leq$ \textsf{Nsteps} ) 
	\begin{enumerate}
	\item Calcular el elemento de la malla que contiene a la
	part\'{\i}cula $\rightarrow$ \textsf{Element}
	\item Calcular la velocidad en \textsf{Position}
	\item Calcular la nueva posici\'on y actualizar \textsf{Position}
	\item \textsf{s = s + 1}
	\end{enumerate}
	{\bf\textsf{END}}
\end{enumerate}
}}

%\emp }

\strut

En el primer paso se determinan las condiciones iniciales para comenzar
el c\'alculo. Las mallas sobre las cuales se puede aplicar el algoritmo
pueden ser de alguno de los tipos mostrado en la figura \ref{figmallas}.
En el paso 1.{\it{d)}} se determina la posici\'on inicial. Debido a
que el campo vectorial no cambia en el tiempo, se lee un s\'olo campo
vectorial, paso 2. Finalmente, en el paso 3, es donde se implementan los 
algoritmos num\'ericos importantes para el c\'alculo de la streamline:
3.{\it{a)}} {\bf{Localizaci\'on}};
3.{\it{b)}} {\bf{Interpolaci\'on}} o evaluaci\'on de una funci\'on;
3.{\it{c)}} {\bf{Integraci\'on}}.
Estos c\'alculos se describen en la secci\'on 3.


\subsection{Pathlines, streaklines y timelines}
Para flujos que dependen del tiempo, se pueden definir {\bf{pathlines}},
{\bf{streaklines}} y {\bf{timelines}}. Las primeras son trayectorias
de part\'{\i}culas liberadas desde posiciones fijas.
Una streaklines se construye conectando mediante una curva y en un instante
de tiempo fijo, las posiciones de un conjunto de part\'{\i}culas que han
sido previamente liberadas, en instantes de tiempo diferentes, desde la
misma posici\'on inicial, v\'ease figura \ref{streaktime}.  
Las timelines son l\'{\i}neas que conectan un conjunto de part\'{\i}culas
soltadas en diferentes posiciones iniciales y en el mismo instante de
tiempo, v\'ease figura \ref{streaktime}.   
En todos los casos, el campo de vectorial var\'{\i}a con el tiempo y la
$\Delta s$, usada para las streamlines, se sustituye por la $\Delta t$,
que es el incremento de tiempo con el que los campos vectoriales son
calculados. Por supuesto, la $\Delta t$ puede ser variable.
Existe entonces, un lapso de tiempo bien definido, durante el cual el campo
vectorial es v\'alido. Por lo tanto, en cada paso de tiempo, se debe leer
o calcular un nuevo campo vectorial.  

En este caso se debe integrar la ecuaci\'on (\ref{ode}), tomando en cuenta
que la velocidad depende del tiempo. 
El algoritmo general para calcular pathlines es como sigue:

\strut

%\framebox{ \bmp{5.5in}

\begin{center}\fbox{\bf{Algoritmo 2: Pathlines}}\end{center}

{\small{
\begin{enumerate}
\item Se determina lo siguiente:
	\begin{enumerate}
	\item Tipo de malla sobre la que est\'a definido el campo
	vectorial $\rightarrow$ \textsf{Mesh}
	\item N\'umero de pasos en el tiempo para calcular la trayectoria 
	$\rightarrow$ \textsf{Nsteps}
	\item Tama\~no del paso en el tiempo ($\Delta t$) para calcular la
	trayectoria $\rightarrow$ \textsf{Dstep}
	\item Posici\'on inicial de la part\'{\i}cula 
	$\rightarrow$ \textsf{Position}
	\end{enumerate}
\item {\bf\textsf{WHILE}} ( \textsf{Position} $\in$ \textsf{Mesh}
			      \quad {\bf\textsf{OR}} \quad 
			    \textsf{t} $\leq$ \textsf{Nsteps} ) 
	\begin{enumerate}
	\item Obtener el campo vectorial $\rightarrow \vec{U}(t)$
	\item Calcular el elemento de la malla que contiene a la
	part\'{\i}cula 
	\item Calcular la velocidad en \textsf{Position}
	\item Calcular la nueva posici\'on (actualizar \textsf{Position})
	\item \textsf{t = t + 1}
	\end{enumerate}
	{\bf\textsf{END}}
\end{enumerate}
}}

%\emp }

\strut

Como se puede observar, hay un cambio m\'{\i}nimo pero importante con respecto 
al algoritmo de streamlines. 
En este caso, el campo vectorial se actualiza en cada paso.
Los algoritmos para calcular streaklines y timelines son muy similares al
del c\'aclulo de pathlines. Para las streaklines, se debe determinar el
intervalo de tiempo para soltar una nueva part\'icula desde la posici\'on
inicial. Adem\'as, para cada part\'{\i}cula, se debe aplicar el ciclo del
paso 2 y unir las part\'{\i}culas en los instantes de tiempo deseados. 
En el caso de las timelines, se deben determinar puntos iniciales para
construir una l\'{\i}nea y el ciclo del paso 2, se aplica a cada punto de la
l\'{\i}nea.  

\begin{figure}[h!]
$$\includegraphics[height=8cm]{streak.eps}$$
\caption{Streaklines y Timelines}
\label{streaktime}
\end{figure}

Una generalizaci\'on de las timelines consiste en definir una superficie
en tres dimensiones y aplicar los algoritmos 1 o 2, seg\'un sea el caso,
para generar superficies deformadas. Con el objeto de ser consistentes,
nombraremos a ese tipo de seguimiento como timesurfaces v\'ease 
figura \ref{surface}. Para obtener m\'as informaci\'on respecto de 
seguimiento de superficies v\'ease \cite{}.

\begin{figure}[h!]
$$\includegraphics[height=5cm]{surface.eps}$$
\caption{Timesurfaces}
\label{surface}
\end{figure}

\section{Algoritmos num\'ericos}

En los algoritmos 1 y 2, hay pasos en donde se realizan c\'alculos
num\'ericos importantes. Se debe elegir el algoritmo correcto para
disminuir el error y mejorar el rendimiento de la implementaci\'on. En
seguida se discuten algunos de los m\'etodos m\'as empleados para estos
c\'alculos. 

\subsection{Localizaci\'on}

La localizaci\'on consiste en determinar la posici\'on de la part\'{\i}cula
dentro del dominio. En el caso en que el dominio est\'e definido mediante
una malla regular estructurada uniforme, este proceso se transforma en 
encontrar el elemento que contiene a la part\'{\i}cula, mediante el siguiente
algoritmo: 

\begin{enumerate}
\item Se dividen los valores de las coordenadas de la posici\'on de la 
part\'{\i}cula, entre los incrementos de la malla en cada una de sus 
direcciones:
\begin{equation}
(a, b, c) = (\frac{x}{\Delta x}, \frac{y}{\Delta y}, \frac{z}{\Delta z})
\end{equation}

\item Se separan la parte entera y decimal de la operaci\'on anterior:

\begin{equation}
(a, b, c)=[i,j,k]+(\alpha ,\beta ,\gamma), 
\end{equation}

\noindent
donde $i,j,k \in I$ y $\alpha ,\beta ,\gamma \in [0,1]$, v\'ease figura 
\ref{celda}.

\item La parte entera $[i,j,k]$ nos dice en que elemento de la malla se 
encuentra la part\'{\i}cula mientras que los decimales 
($\alpha ,\beta ,\gamma$), se utilizan para la interpolaci\'on.
\end{enumerate}


\begin{figure}[h!]
$$\includegraphics[height=6cm]{celda.eps}$$
\caption{Partes entera y fraccional de las coordenadas de la posici\'on de
la part\'{\i}cula.}
\label{celda}
\end{figure}

\subsection{Interpolaci\'on}

Si $\vec{x}$ esta dentro del elemento de malla $[i,j,k]$ y tiene las partes
fraccionales $(\alpha ,\beta ,\gamma)$, entonces, una f\'ormula de
interpolaci\'on para el campo vectorial es como sigue: 


{{ 
 \begin{eqnarray} 
 \vec{U}(i+\alpha, j+\beta, k+\gamma) & = &
 \Big{[}\Big{(}\vec{U}_{i,j,k}(1-\alpha)+\vec{U}_{i+1,j,k}\alpha\Big{)}
 \Big{(}1-\beta\Big{)} + 
 \nonumber \\ 
& & \Big{(}\vec{U}_{i,j+1,k}(1-\alpha)+\vec{U}_{i+1,j+1,k}\alpha\Big{)} \beta
 \Big{]}\Big{(}1-\gamma\Big{)}+ \nonumber \\
& &
 \Big{[}\Big{(}\vec{U}_{i,j,k+1}(1-\alpha)+\vec{U}_{i+1,j,k+1}\alpha\Big{)}\Big{(}1-\beta\Big{)}+  
 \nonumber \\
& &
 \Big{(}\vec{U}_{i,j+1,k+1}(1-\alpha)+\vec{U}_{i+1,j+1,k+1}\alpha\Big{)}\beta\Big{]}\gamma 
 \end{eqnarray}
}}

 \noindent
 donde $\vec{U}_{i,jk}$ es el valor del campo vectorial definido en el
 punto de la malla $(i,j,k)$, v\'ease figura
 \ref{celda}. Esta f\'ormula se conoce como interpolaci\'on trilineal y es
 r\'apida y simple, pero la exactitud puede perderse cuando 
 la malla es gruesa o el dominio es irregular. 
 En Yeung and Pope \cite{Yeung} se comparan f\'ormulas de interpolaci\'on 
 de orden m\'as alto. Por ahora se utilizar\'a la interpolaci\'on trilineal
 para todos los c\'alculos.

\subsection{Integraci\'on}

La integraci\'on de la ecuaci\'on (\ref{ode}) es muy impotante, pues 
es aqu\'{\i} donde se obtienen los puntos que componen la trayecttoria. Los 
errores num\'ericos mas dr\'asticos pueden ser introducidos en esta parte del
c\'aculo. La integraci\'on se realiza de la siguiente manera: sea
$\vec{x}(t)$ la posici\'on inicial de la part\'{\i}cula, entonces,
despu\'es de un cierto $\Delta t$ la nueva posici\'on ser\'a:

\begin{equation}\label{integra}
\vec{x}(t + \Delta t) = \vec{x}(t) + \int_t^{t + \Delta t} \vec{U}(\vec{x},t) dt
\end{equation}

El segundo t\'ermino de la derecha de la ecuaci\'on (\ref{integra})
determina el m\'etodo de integraci\'on y el error introducido. 

\subsubsection{Euler}

El m\'etodo de Euler es el m\'as simple de todos. Se considera que la
velocidad se mantiene constante durante el intervalo $[t, t+\Delta t]$,
entonces la ecuaci\'on (\ref{integra}) se transforma en:

\begin{equation}\label{Eintegra}
\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{U}(\vec{x},t) \Delta t
\end{equation}

\noindent
donde del lado derecho de la ecuaci\'on se tienen los primeros dos
t\'erminos de una serie de Taylor. Obviamente, este m\'etodo tiene un error
del orden de $\Delta t$. Cuando $\Delta t$ es peque\~no el error puede no
ser significativo. Sin embargo, no tenemos control sobre $\Delta t$ y es
posible que el error sea muy grande.

\subsubsection{Runge-Kutta de segundo orden (mid-point)}

Este m\'etodo (tambi\'en conocido como m\'etodo de Heun) consiste de los
siguientes pasos.

\strut

Sea $\vec{x}(t) = \vec{x}_k$ y $\vec{x}(t+\Delta t) = \vec{x}_{k+1}$, 
entonces:

 \begin{eqnarray*}
 a & = & \vec{U}(\vec{x}_k, t) \Delta t, \\ \nonumber 
 b & = & \vec{U}(\vec{x}_k + \frac{a}{2}, t+\frac{\Delta t}{2}) \Delta t, \\ \nonumber 
 \vec{x}_{k+1} & = & \vec{x}_{k}+ b, \\ \nonumber
 t & = & t + \Delta t, \\ \nonumber
 k & = & k+1, \\ \nonumber
 \end{eqnarray*}

El grado de aproximaci\'on es del orden de $\Delta t ^2$.

\subsubsection{Runge-Kutta de cuarto orden}

Un m\'etodo m\'as exacto es conocido como Runge-Kutta de cuarto orden. El
m\'etodo consiste de los siguientes pasos:

\noindent
 Sea $t = t_{1}$, $k = 0$ y $\vec{x}_k$ la posici\'on actual de la
 part\'{\i}cula. Entonces, para $t < t_{n}$ y $\vec{x}_{k}$ dentro de la
 malla hacer lo siguiente:
 
 \begin{eqnarray*}
 a & = & \vec{U}(\vec{x}_{k},t) \Delta t, \\ \nonumber
 b & = & \vec{U}(\vec{x}_{k}+a/2, t + \frac{\Delta t}{2}) \Delta t,\\\nonumber
 c & = & \vec{U}(\vec{x}_{k}+b/2, t + \frac{\Delta t}{2}) \Delta t,\\\nonumber
 d & = & \vec{U}(\vec{x}_{k}+c) \Delta t + \Delta t, \\ \nonumber
 \vec{x}_{k+1} & = & \vec{x}_{k}+(a+2b+2c+d)/6, \\ \nonumber
 t & = & t + \Delta t, \\ \nonumber
 k & = & k+1, \\ \nonumber
 \end{eqnarray*}

El grado de aproximaci\'on es del orden $\Delta t ^4$,
por lo que el error resulta bastante peque\~no.


\section{Emisores}

La posici\'on inicial de una trayectoria, que es la condici\'on inicial
para resolver la ecuaci\'on (\ref{ode}), tambi\'en se conoce como punto
emisor o simplemente emisor. En algunos casos se tiene un conjunto de
emisores que pueden o no describir una forma geom\'etrica. Los emisores
pueden ser fijos o m\'oviles. En el \'ultimo caso se debe tener la
facilidad de modificar la posici\'on del emisor de manera interactiva (usando
el teclado, mouse, wand, etc.) y precisa.

\subsection{Emisores simples}

En este caso los emisores se calculan de una manera simple y pueden ser:

\begin{itemize}
\item Un punto.
\item Un arreglo de puntos definidos en una malla regular en 1D (l\'{\i}nea), 
2D (plano) o 3D (cubo).
\item Un arreglo de puntos definidos sobre formas geom\'etricas (c\'{\i}rculos,
esferas, elipses, elipsoides, pol\'{\i}gonos, etc.).
\end{itemize}

En estos casos s\'olamente es necesario conocer la posici\'on de los puntos.

\subsection{Emisores conexos}

Un emisor conexo consiste de un arreglo de puntos, definido como se mencion\'o
en la secci\'on anterior, y de su conectividad. La conectividad indica los
puntos que se conectan para formar elementos, los cuales pueden ser 
triangulares o cuadrangulares. La figura \ref{conectividad} muestra un 
ejemplo de conectividad.

\begin{figure}[h!]
$$\includegraphics[height=4cm]{conectividad.eps}$$
\caption{Para definir los elementos se utiliza la conectividad. 
(a) Conectividad para cuadr\'angulos: (0,1,5,4), (1,2,6,5),
(2,3,7,6), etc. (b) Conectividad para tri\'angulos: (0,1,5), (0,5,4), (1,2,6),
(1,6,5), (2,1,6), (1,7,6), etc. 
Obs\'ervese que siempre se sigue la misma orientaci\'on (contraria a las 
manecillas del reloj).}
\label{conectividad}
\end{figure}

En este caso la conectividad es importante, pues se requiere de reconstruir una
superficie en tiempos subsecuentes que muestre deformaciones de acuerdo con
el campo vectorial que se est\'e analizando, v\'ease figura \ref{surface}.

\section{An\'alisis y Dise\~no (draft)}

En esta secci\'on se describe el dise\~no y la forma de implementaci\'on del
sistema para Seguimiento Lagrangiano de Part\'{\i}culas. La idea es que este
sistema sea modular y adaptable a diferentes aplicaciones. Se busca baja
cohesi\'on y alto nivel de abstracci\'on. El an\'alisis, dise\~no y 
notaci\'on se basan en el modelo UML, v\'ease \cite{}.


\subsection{Glosario de clases}

El primer glosario de clases, de acuerdo con los requerimientos
\footnote{Los requerimientos en este caso se obtienen de las secciones
2,3 y 4, en donde se describe la teor\'{\i}a y componentes de este
sistema. Estos se modificaran conforme avanze el proyecto.}
que se definen en este proyecto es el siguiente:
\footnote{Los nombres de las clases estar\'an en ingl\'es.}

\begin{description}
\item[Trace]: Representa una trayectoria de una part\'{\i}cula en general.
  Los tipos particulares de trayectorias que se pueden calcular (y que
  en principio heredan de Trace) son:
  \begin{description}
  \item[StreamLine]: Definida en el texto.
  \item[PathLine]: Definida en el texto.
  \item[StreakLine]: Definida en el texto.
  \item[TimeLine]: Definida en el texto.
  \item[TimeSurface]: Definida en el texto.
  \end{description}

\item[Particle]: Representa a una part\'{\i}cula.
\item[VectorField]: Representa el campo vectorial.
\item[Domain]: Representa el dominio de c\'alculo.
\item[Interpolation]: Para interpolar la velocidad. Existen diferentes
 maneras de hacer una interpolaci\'on:
  \begin{description}
  \item[Trilineal]: Definida en el texto.
  \item[Cuadr\'atica]: V\'ease \cite{}.
  \item[B-splines]: V\'ease \cite{}.
  \end{description}
\item[Solver]: Contiene m\'etodos para integrar la ecuaci\'on (\ref{ode}).
  Existen diferentes m\'etodos de interpolaci\'on:
  \begin{description}
  \item[Euler]: Definido en el texto.
  \item[Runge-Kutta $2^o$ orden]: Definido en el texto.
  \item[Runge-Kutta $4^o$ orden]: Definido en el texto.
  \end{description}

\item[Emitter]: Representa a los diferentes tipos de emisores.
  Los diferentes tipos de emisor son:
  \begin{description}
  \item[Simple]: Definido en el texto.
  \item[Connected]: Definido en el texto.
  \end{description}
\end{description}


\subsection{Diagrama de clases}

\begin{figure}[h!]
$$\includegraphics[height=3cm]{traces.eps}$$
\caption{}
\label{traces}
\end{figure}

\begin{figure}[h!]
$$\includegraphics[height=2cm]{interpolacion.eps}$$
\caption{}
\label{interpolacion}
\end{figure}

\begin{figure}[h!]
$$\includegraphics[height=2cm]{solver.eps}$$
\caption{}
\label{solver}
\end{figure}

\begin{figure}[h!]
$$\includegraphics[height=2cm]{emisor.eps}$$
\caption{}
\label{emisor}
\end{figure}

%\begin{thebibliography}{99}


 %\bibitem{distresf} ...

 %\bibitem{Hall} M. Hall and J. Warren, Adaptive Polygonization of
 %Implicitly Defined Surfaces, IEEE Computer Graphics and Applications,
 %{\bf{10}}, 33-42, 1990.

%\bibitem {uml}
%Jacobson, I., Booch, G. y Rumbaugh, J.,
%{\it{The Unified Software Development Process}},
%Addison-Wesley,
%{\bf{1999}}.

% \bibitem{Yeung} P. Yeung and S. Pope, An Algorithm for Tracking Fluid
% Particles in Numerical Simulations of Homogeneous Turbulence,
% \emph{J. of Comp. Physics}, {\bf{79}}, 373-416, 1988.

%\end{thebibliography}


\end{document}















